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We present a wide array of quantum measures on numerical solutions of ID Bose- and Fermi- 
Hubbard Hamiltonians for finite-size systems with open boundary conditions. Finite size effects 
are highly relevant to ultracold quantum gases in optical lattices, where an external trap creates 
smaller effective regions in the form of the celebrated "wedding cake" structure and the local density 
approximation is often not applicable. Specifically, for the Bose-Hubbard Hamiltonian we calculate 
number, quantum depletion, local von-Neumann entropy, generalized entanglement or Q-measure, 
fidelity, and fidelity susceptibility; for the Fermi-Hubbard Hamiltonian we also calculate the pairing 
correlations, magnetization, charge- density correlations, and antiferromagnetic structure factor. Our 
numerical method is imaginary time propagation via time-evolving block decimation. As part of our 
study we provide a careful comparison of canonical vs. grand canonical ensembles and Gutzwiller 
vs. entangled simulations. The most striking effect of finite size occurs for bosons: we observe a 
strong blurring of the tips of the Mott lobes accompanied by higher depletion, and show how the 
location of the first Mott lobe tip approaches the thermodynamic value as a function of system size. 



INTRODUCTION 



Quantum phase transitions (QPTs) are often treated 
in the thermodynamic limit, where the concepts of a crit- 
ical point, critical exponents, and scaling relations are 
very clear However, many physical many-body sys- 
tems for practical applications these days are finite. For 
example, nano-devices can be smaller than 100 nm. Then 
a crystal with a lattice constant of a fraction of a nanome- 
ter has just 2 or 3 orders of magnitude to work with 
in scale. In nuclear physics one has QPTs in finite nu- 
clei which are very far from being in the thermodynamic 
limit Experiments on trapped ions, recently used to 
simulate many body Hamiltonians, normally work with 
small systems [3] . Ultracold atoms in optical lattices pro- 
vide a fourth example, where the number of lattice sites 
is typically 50 to 100 in each direction, so that in one di- 
mension (ID) they are highly mesoscopic. These systems 
have been proposed as quantum simulators to solve hard 
quantum many body problems, including the positive- 
Fermi-Hubbard Hamiltonian for general filling [3, [H, @, 0] • 
Progress towards this goal requires understanding of how 
finite-size effects change the quantum phase diagram. In 
this paper, we explore such mesoscopic or finite size ef- 
fects for both Bose- and Fermi-Hubbard Hamiltonians. 

In addition to finite-size effects, ultracold atomic quan- 
tum simulators are non-uniform, as ordinarily there is a 
weak overall harmonic trap keeping atoms from spilling 
off the edge of the lattice. This non-uniformity gives rise 
to a "wedding-cake" or tiered structure of Mott plateaus 
separated by superfluid layers for bosons l^l^, and sim- 
ilar layered structure for fermions |ll|, [l2|. The non- 
uniformity creates even smaller regions of Mott insulator 
or superfluid, often just 10 or 20 sites in length. For 
these small regions the quantum phase diagram is sub- 
stantially different than in the thermodynamic limit: for 
example, the quantum depletion, (defined precisely be- 
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FIG. 1: (color online) The first Mott lobe of the Bose-Hubbard 
Hamiltonian. Quantum depletion for a system of interact- 
ing bosons on 51 lattice sites at T = as a function of 
chemical potential /x and hopping parameter t, both normal- 
ized to interaction strength [/, obtained with Time-Evolving 
Block Decimation (TEBD) in imaginary time propagation, as 
discussed below. The white circular points represent best- 
converged results of Density Matrix Renormalization Group 
methods for the boundary of the Mott lobe in the thermody- 
namic limit Nl. 



low in Sec.[TT|) shows that the Mott lobe is not closed but 
instead goes through a very broad crossover into the su- 
perfluid regime, as we show in Sec. lIII Al Figure [T] shows 
the sharp phase boundary of a Bose Hubbard system in 
the thermodynamic limit (white circles Q) and the con- 
tinuously variable quantum depletion (defined in Sec. HH) 
obtained in our calculations for 51 sites. As can be seen 
in the figure, the tip of the Mott lobe is significantly 
blurred even for this relatively large lattice. 
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Approximating the effects of a slowly varying trap 
by using a local chemical potential in the local density 
approximation, fi = /i(r), is not very useful for cold 
atoms in ID. The breakdown of the local density approx- 
imation has been explored with quantum Monte Carlo 
and exact diagonalization elsewhere (see [13] and refer- 
ences therein); our work, which isolates finite size effects 
in small uniform regions, is complementary to previous 
studies which have explored the harmonic trap structure 
as a whole [ol, [isl, [l3 • The success of quantum simulators 
relies in part on understanding all possible sources of er- 
ror in the same sense as one might do so for a high 
precision measurement of time-reversal symmetry, for ex- 
ample. Therefore going beyond the local density approx- 
imation is key to the eventual success of quantum 
simulators, especially in lower dimensions. Our use of 
TEBD offers a wide array of quantum measures, many 
of which are not directly accessible to quantum Monte 
Carlo [3 or exact ID methods [19, 20], includ- 

ing various entanglement entropies and fidelity. Multiple 
measures provide maximal information, useful for a thor- 
ough analysis; in particular, fidelity and fidelity suscepti- 
bility make no assumptions about states or symmetries, 
but simply ask about many body wavefunction overlap 
in a model-independent way [2l| . 

It has already been shown that the mean-field 
Gutzwiller approximation, which assumes a product 
state between sites of the lattice, is not very useful in 
ID for Hubbard Hamiltonians [6]. But how far does one 
have to move away from the Gutzwiller approximation 
in order to get correct results? A high level of effort 
has been put into alternative methods, such as gener- 
alized dynamical mean field theory in two dimensions 
for Bose- Fermi mixtures [22]. The convergence param- 
eter in TEBD, which can be couched in terms of the 
degree of spatial entanglement, allows one to go beyond 
the Gutzwiller approximation in a methodical manner, 
and we present a thorough study of this issue. A second 
question concerning numerical methods is the use of the 
canonical ensemble vs. the grand canonical ensemble. 
The two only formally give the same result in the ther- 
modynamic limit. The canonical ensemble is numerically 
preferred for TEBD and related methods (2^ because 
the number of Fock states is much smaller and there- 
fore calculations are more efficient. We explicitly explore 
the effects of ensemble dependence Hubbard Hamiltonian 
phase diagrams. 

This article is outlined as follows. In Sec.[lTlwe present 
the Bose- and Fermi-Hubbard Hamiltonians, the many 
quantum measures we use, and a brief description of 
and motivation for our implementation of TEBD. In 
Sec. mil we present and discuss our numerical results 
for mesoscopic effects in ground states of the Bose- and 
Fermi-Hubbard Hamiltonians, respectively. In Sec. IIVI 
we compare canonical to grand canonical ensembles and 
Gutzwiller to entangled TEBD simulations; these com- 
parisons also serve as a convergence study. Finally, in 
Sec. we conclude. 



II. MODELS, MEASURES, AND METHODS 

We treat two fundamental Hamiltonians for ultracold 
quantum gases in optical lattices, the Bose- and Fermi- 
Hubbard Hamiltonians. In both cases we take the tight- 
binding lowest-band approximation. The Bose-Hubbard 
Hamiltonian is 

i 

where t is the hopping strength, U is the interaction 
strength, /i is the chemical potential, 6^, h\ are bosonic 
destruction/creation operators satisfying bosonic com- 
mutation relations on-site and commuting between sites, 
n^^^ is the bosonic number operator, and (i, j) denotes a 
sum over nearest neighbors. The corresponding Fermi- 
Hubbard Hamiltonian in particle-hole symmetric form is 

i 

where fi^ , //^ are fermionic destruction and creation op- 
erators satisfying fermionic anticommutation relations; t 
and U have the same meaning as in the bosonic case; and 
we have assumed s-wave interactions and spin 1/2, so 
that (J G {T7 i}- We note that there is an exact mapping 
under the canonical transformation from the negative- 
to positive- 1/ Fermi-Hubbard Hamiltonian at half filling 
(/i = 0). This will appear visually in our phase diagrams 
in Sec. HITBl 

To find the ground states of these models, we use 
TEBD [23] in ID and propagate in imaginary time, start- 
ing with a state with uniform weights. TEBD has one 
main convergence parameter, the number of eigenvalues 
X retained in the reduced density matrix. In all cases 
we have checked convergence by at least doubling x, and 
seeing if we obtain the same results in the plots; all plots 
have at least a x of 16, and we have spot-checked a x of 
32. We emphasize that as we have a very high resolu- 
tion of (250 X 100) pixels in our diagrams, we cannot go 
to the preferred typical values of x in the hundreds, in 
common use for TEBD calculations; also, grand canon- 
ical calculations are inherently more restrictive, due to 
the much larger state space. We use both canonical and 
grand canonical ensembles depending on circumstances. 
The latter is useful to understand the proper interpo- 
lation method to be used for the canonical ensemble 
for smaller numbers of sites. In the canonical ensem- 
ble, chemical potential is derived by finite differencing, 
/i ^ E{N + 1) — E{N); for just a few sites, this would 
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lead to only a few points of resolution on the chemical 
potential axis without interpolation. That is, a finite dif- 
ference approximation for the chemical potential is only 
directly useful for large TV; interpolation can only be jus- 
tified by comparing to the grand canonical results. We 
illustrate this point in Sec. IIV Ai For fermions, there is 
an additional advantage of grand canonical calculations 
in that, in practice, we find that imaginary time propa- 
gation is less likely to get stuck in a local minimum. 

For fermions, the local, or on-site dimension of the 
Hilbert space is fixed at d = 4 since we consider only 
a single band. For bosons d provides a second conver- 
gence parameter. In practice we truncate at d = 5 + 1, 
meaning to 5 bosons per site, and we have checked 
(i = 7 + 1 for specific points in the phase diagrams and 
found no visible changes. We implement imaginary time 
propagation in the second order Suzuki- Trotter scheme 
to obtain the ground state. Finally, we parallelize over 
the parameter space (t/U^ji/U) via MPI for bosons, and 
{U/t^li/t) for fermions. We typically perform our calcu- 
lations on between 5x8 and 60 x 8 cores. We note that 
TEBD becomes exact, aside from the O {pt^) error of the 
Suzuki- Trotter decomposition, for x ^ d^^^'^K 

Our quantum measures can be divided into four 
classes: moments, correlation, entropies, and fidelity. Al- 
though these classes of measures can have significant 
overlap in the information they provide, nevertheless 
there are always cases where they are distinct. For ex- 
ample, the spin singlet is maximally entangled but is not 
correlated. Thus it is useful to consider all classes. 

First, for moments, the energy, £^5, and filling factor, 
are defined as 



The charge-density wave correlation, C, is given by 

c = {ee^), (8) 

(9) 
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The antiferromagnetic structure factor is given by 



5f 



(10) 



(11) 



Third, for entropies, the average local von Neumann 
entropy [1^, or average local entropy of entanglement, 
S, is 



1 ^ 

S = -iJ2^' (pilogrfPi) G [0,1] 



(12) 



where the local reduced density matrix, p^, is defined as 



Pi 



Tr 



(13) 



and p is a pure-state density matrix for the whole system, 
as we work at zero temperature. Note that the base 
d of the log normalizes the entropy to lie between zero 
and unity. We also studied the Q-measure [1^, [sO] or 
average local impurity, Q, a special case of generalized 
entanglement [3l|, [32[ , which gives similar results to the 
von Neumann entropy (see Sec. IIIip : 



(3) 
(4) 



where L is the number of lattice sites. To calculate Ef 
and Vf for fermions, one must also sum over the spin 
index a. 

Second, for correlations in bosons we calculated the 
quantum depletion, which is a measure of how far 
away the system is from the Penrose- Onsager definition 
of superfiuidity [25[: 



D = 1 



Ai 



(5) 



where {A^} are the eigenvalues of the single particle den- 
sity matrix {blbj) ordered so that Ai > A2 > • • • > A^. 
For correlations in fermions we use three standard mea- 
sures [13 . The s-wave pairing correlation, P, is given 

by 



P = (AtA + AAt), 

At.^E/t/^ 



(6) 
(7) 



d 



d-1 



l-iETr(P?) 



e [0, 1] . (14) 



Fourth, we calculated both fidelity, /, and fidelity sus- 
ceptibility [21I, [3I, [33, Xf- The fidelity is an overlap 
measure on the states, and is given by 



/(C,C + <^C)^ 1(^(0 l^(C + <^C))l , 



(15) 



where ^ is a control parameter in a quantum phase dia- 
gram, in our case either chemical potential or hopping. 
The fidelity susceptibility is a second derivative of the 
fidelity in the control parameter, and we normalize it to 
the unit length: 



(16) 



This proves to be a more useful quantity than the fidelity, 
as it is independent of the control parameter step size 
and a divergence of Xf the thermodynamic limit 
signals a quantum phase transition. Note that Xf bears 
no relation to the convergence parameter for TEBD, x- 
In the figures below we take ( to be t/U for the Bose- 
Hubbard Hamiltonian and U/t for the Fermi-Hubbard 
Hamiltonian. 
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III. QUANTUM MEASURES ON GROUND 
STATES OF HUBBARD HAMILTONIANS 

A. Bose-Hubbard Hamiltonian 

The Bose-Hubbard Hamiltonian phase diagram has 
been calculated in the mean field approximation in higher 
dimensions [H and with both the DMRG method g [s^ 
and Monte Carlo in one dimension. Trapped Bose sys- 
tems with an external harmonic [13] and linear poten- 
tial [s^l have been treated with Monte Carlo and TEBD, 
respectively. However, these studies did not focus on fi- 
nite size effects without aiming at the thermodynamic 
limit. 

In Fig.[2]are illustrated a suite of quantum measures for 
L = 3, 5, 10, 20, and 40 sites, from left to right in column 
form. The rows represent the different quantum mea- 
sures as labeled on the right hand side: from top to bot- 
tom these are density, entropy, Q-measure, depletion, fi- 
delity, and fidelity susceptibility. The calculation of these 
measures and their definitions are given in Sec. [Til We 
emphasize that quantum Monte Carlo cannot be used to 
obtain entanglement. Only matrix-product-state-based 
techniques can do so for general local Hamiltonians, in- 
cluding TEBD and DMRG [s^ , which lead to an equiva- 
lent result, although with different numerical implemen- 
tations. Bethe-ansatz-based techniques have been used 
to calculate entanglement in certain limitin g ca ses (near 

= 0) in the Fermi-Hubbard Hamiltonian [20], but not 
generally. The case of L = 3 is exact diagonalization for 
X = 15 and d = 6, but at L = 5 the result is already 
approximate, since 6^ = 36 states would be required for 
exact diagonalization in matrix product state form (see 
Sec. [ni). Each panel of Fig. [2] is of resolution 80 x 240 in 
t/U X (i/U. The largest system of L = 40 took several 
weeks to complete in the grand canonical ensemble on 
20 nodes with separate memory, each node consisting of 
8 cores with shared memory, all on the Golden Energy 
Computing Organization [39] high performance comput- 
ing cluster "Ra" . For this reason we do not go to x = 32 
for the whole phase diagram, but only make spot checks, 
since simulations scale as x^- We calculated both canon- 
ical and grand canonical results, and show the canonical 
ones in Fig. [2l 

Turning to the physical interpretation of Fig. [2l we 
begin with number. In the canonical ensemble a definite 
total number is guaranteed. This can be seen first in 
the average filling, which proceeds in multiples of 1/L, 
as illustrated in the top row of the figure, starting from 
the vacuum and a filling of zero. The vacuum region is 
evident in the lower portion of the panels in the top row 
of where it is shaded a deep blue. The density increases 
in discrete jumps as the chemical potential is increased; 
interpolation is used between curves of definite number 
difference in order to fill in the picture, as described in 
Sec. [HI and Sec. [IVA[ 

In the second row we show the quantum depletion. One 
clearly sees the highly depleted Mott lobes and the lowly 



depleted outer superfiuid region. Note that there are only 
two Mott lobes displayed in these figures; the apparent 
third, lowest such lobe is in fact the vacuum, = 0, for 
which quantum depletion is ill-defined. In keeping with 
the ID nature of the problem we consider, we observe 
that even the superfiuid region is still depleted at a level 
of about 30 % far away from the Mott lobes. Again, the 
smooth nature of the Mott-insulator-superfiuid transition 
near the tip is in strong contrast to the sharper transition 
in other regions. 

In the third and fourth rows of Fig. [2] we show en- 
tropies. For completeness we compare the von Neu- 
mann entropy or entropy of entanglement, which is now 
a standard measure used in quantum information [13] 
and also proposed for the study of quantum phase tran- 
sitions [2^, [41I, [42[, to the Q-measure, which has been 
more closely connected to local observables [3l|, |43| . We 
find both measures useful and nearly but not exactly the 
same; we choose the Q-measure for the rest of our study. 
Note that the local von Neumann entropy of the Mott 
insulating regions is low by this measure, which refiects 
the nearly pure number state on each site in the Mott 
phase. On the other hand, the superfiuid region displays 
high entropy, because the superfiuid state is delocalized 
across sites, leading to mixed states on each site when 
measured separately. This contrasts with the usual asso- 
ciation of low thermal entropy with a superfiuid phase. 

In the second through fourth rows one can clearly pick 
out the Mott-insulating region. It has a sharp phase 
boundary for small t/U < 0.2. Kiihner, White and 
Monien [8] first provided an accurate evaluation of the 
critical value, tc = 0.29 ±0.1, for the Kosterlitz-Thouless 
(KT) transition at the tip of the Mott lobe in the infi- 
nite system limit; in our notation this corresponds to a 
value of t/U = 0.29. Near tc the entropy changes more 
continuously for 3, 5, and 10 sites, while for 20 sites the 
lobe begins to close. Finally at 40 sites, on the far right 
of the figure, the region of continuous change is very nar- 
row. We clarify that by "sharper boundary" as compared 
to "more continuous" we are speaking in relative terms, 
specifically of the difference in entropies between states 
with commensurate and incommensurate filling. There 
is no truly sharp transition since this is a finite system; 
all "quantum phase transitions" seen here are in fact for- 
mally crossovers between regions with different charac- 
teristics. 

Finally, in the fifth and sixth rows we show the fidelity 
/ and fidelity susceptibility Xf- present imple- 

mentation we compare the neighboring pixels of the fig- 
ures in t/U for curves of fixed N and L, as described 
in Sec. [Ill This excludes a comparison between changes 
in total number; such a comparison, which we also per- 
formed in the grand canonical ensemble, is not shown, 
because features in / and x / are dominated by changes in 
the conserved quantity. In the thermodynamic limit the 
superfiuid region of the phase diagram forms a continuous 
tensor product space and the fidelity is everywhere zero 
even for arbitrarily small S(. This is the well-known An- 
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FIG. 2: (color online) Quantum Measures for the Bose-Hubbard Hamiltonian. Shown is the density (first row), average local 
entropy of entanglement (second row), generalized entanglement or Q-measure (third row), quantum depletion (fourth row), 
fidelity (fifth row), and fidelity susceptibility (sixth row). The panels show increasing numbers of sites from left to right: 3, 5, 
10, 20, and 40. 



derson orthogonality catastrophe [i^. Thus, the fidelity, 
through its scaling and qualitative behavior, reflects for 
a finite system the critical behavior of an infinite sys- 
tem. Indeed, the Mott lobes appear clearly in this set of 
measures, as can be seen in the figure. 

In order to provide a more quantitative estimate of the 



tip location for finite systems than a by-eye estimation 
from Fig. [21 we study quantum measures and derivatives 
thereof for unit filling = L in the canonical ensemble. 
Three figures illustrate our results. In Fig. [3] are shown 
derivatives of quantum measures which could be used to 
define a crossover; in rows from top to bottom these are 
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derivatives of average local von Neumann entropy, Q- 
measure, and quantum depletion; and lastly, the fidelity 
susceptibility. The fidelity susceptibility provides a par- 
ticularly sharp signature of the Mott phase, a result that 
was noted previously in the context of perodic bound- 
ary conditions and the completely connected graph [33^ . 
Within the n*^ Mott lobe, Xf assumes values close to 
those given by second order perturbation theory [s^] in 
the limit that t/U ^ 0^: 



X/ = 2n(n + l)- 



1 



L 



(17) 



for open boundary conditions. For periodic boundary 
conditions the system size cancels and the expression is 
Xf = 2n(n + 1). These perturbations represent particle- 
hole fiuct nations, which are weighted more heavily in 
higher Mott lobes due to Bose stimulation. This causes 
X/ to have a higher value in the upper Mott lobe in Fig. [2] 
as compared to the lower lobe. 

As the number of sites increases from 3 to 40 the ex- 
trema move towards the thermodynamic limit tc Fig- 
ured! summarizes the most useful extrema of these curves 
in an easy to read single plot showing that quantum 
depletion and fidelity susceptibility appear to converge 
most rapidly to their asymptotic values in the thermo- 
dynamic limit. However, X/ is affected by the boundary 
conditions, as can be seen in Eq. [T71 as well as in Fig. [5l 
where we compare periodic and open boundary condi- 
tions using both exact diagonalization and TEBD; this 
was also a significant issue in DMRG and perturbation 
theory studies, where earlier calculations [3^ obtained 
a less accurate estimate of the critical point due to the 
choice of boundary conditions [8|] . This figure also serves 
as a check on our code. We emphasize that open bound- 
ary conditions are closer to the conditions of ultracold 
atom experiments at the present time. 
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FIG. 3: (color online) Crossover estimators. First derivative 
(left column) and second derivative (right column) with re- 
spect to t/U of various quantum measures for unit filling in 
the canonical ensemble. The vertical line corresponds to the 
known value of tc- 



B. Fermi-Hubbard Hamiltonian 



The Fermi-Hubbard Hamiltonian in ID has an exact 
solution [19, 45] via the Bethe ansatz; in principle, one 
can extract any information from such a solution, includ- 
ing all measures of interest, and the Bethe ansatz has 
been applied to finite as well as infinite domains. Be- 
cause of their computational complexity, Bethe ansatz 
results are not widely used for performing partial traces 
of the type that we consider here [lO]. Therefore it is 
useful to consider finite size effects for a suite of mea- 
sures simultaneously, including measures not previously 
accessed through exact solution techniques. TEBD or 
DMRG methods also enable one to "turn on" non-mean- 
field effects, i.e. entanglement, in discrete steps in the 
form of the parameter x, as we show in Sec. IIVBI 

Figure [6] displays ground state solutions of the Fermi- 
Hubbard Hamiltonian [Eq. ([2])] in the grand canonical 
ensemble for 3, 5, 9, 10, and 20 sites from left to right in 
columns, following the same pattern as Fig. [2] but with 



measures appropriate to fermions. Since we take x = 16, 
for 3 and 5 sites TEBD is equivalent to exact diagonal- 
ization, as d = 4; in fact, we explicitly checked our results 
with a separate exact diagonalization routine. The first 
row shows the density. For positive U we observe that 
all number states are allowed; as in Sec. IHI A[ the grand 
canonical ensemble settles on a total-number-conserving 
state for these small systems due to the discreteness of 
the spectrum. For example, for L = 3 there are seven col- 
ored regions for increasing (i/t and U/t greater than 0, 
of which all regions are visible up to [//t 2; these cor- 
respond to 0,1,. . .,6 fermions. For L = 5 there are eleven 
colored regions, of which 2 are barely visible due to our 
choice of scale, and for higher L we observe the same pat- 
tern. For odd numbers of sites the exact mapping from 
positive U to negative U in Eq. ([2]) at zero chemical po- 
tential does not hold, and therefore there is no symmetry 
around the vertical axis at U/t = for L = 3,5,9, while 
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FIG. 4: (color online) Estimates of the crossover. Extrema of 
the first derivative of the depletion (red squares), the fidelity 
susceptibility (blue triangles), and the second derivative of the 
Q- measure for unit filling in the canonical ensemble (green cir- 
cles). All derivatives are with respect to t/U. The horizontal 
line again corresponds to the limit tc. The curves are a guide 
to the eye. 



for L = 10, 20 we observe this symmetry. For negative 
[/, particles enter only in pairs, excepting at /i = for L 
even, as can be seen from the figure. 

In the second row of Fig. [6] we show the entanglement 
entropy. One clearly observes the particle-hole symme- 
try between negative and positive chemical potential. For 
odd numbers of sites, the most highly entangled regions 
for positive U are those in which there is one extra par- 
ticle or hole above half filling; for even systems it is two 
extra particles or holes, as can be seen in the stripe-like 
pattern for L = 20, for example. States nearer the band 
insulator, with just a few particles (holes) compared to 
the vacuum (maximal filling) are less entangled. Entan- 
glement also decreases in the Mott region for increasing U 
since the system moves closer to a product state. There 
is a discontinuity in the Q-measure between regions of 
half filling and slightly away from half filling which van- 
ishes at [7=0 . This can be understood as the requirement 
that the derivatives approaching from above and below 
are required to satisfy dQ / dn\n^i+ = —dQ/dn\r^=i- by 
particle-hole symmetry. For negative U the most highly 
entangled regions are near zero chemical potential and 
near half-filling. Therefore, in general the pattern for 
both positive and negative U is that the entanglement 
is largest near half filling. We show only the Q-measure, 
having established in Sec. lIII Al that the average local von 
Neumann entropy does not show substantial additional 
information. 

In the third row of Fig. [6] we show the charge-density 
correlation (CDC). This measure shows how close the 
system is to a "checkerboard" occupying every other site. 
For an odd number of sites, as visible in L = 3, 5, 9, and 



FIG. 5: (color online) Open vs. Periodic Boundary Condi- 
tions. Fidelity susceptibility for fixed = L in a compara- 
tive study of open (OBC) and periodic boundary conditions 
(PBC); the latter can only be obtained via exact diagonaliza- 
tion (ED) in our present code, while OBCs can be compared 
directly with TEBD. We note that (1) PBC have a substan- 
tially different fidelity for smaller systems, with a peak at 
much smaller t/U, and (2) TEBD and exact diagonalization 
match up to the O {St^) error from the Suzuki- Trotter de- 
composition for small systems, validating our code. 



negative /7, the CDC is higher for two pairs than for 
one pair, by the nature of the measure [see Eq. (|9])] - 
this is due to the fact that there is one extra pair which 
cannot be placed. For even numbers of sites, as seen 
in L = 10, 20 it is particle-hole symmetric, as expected, 
yielding the same result upon flipping all plots around the 
/i = horizontal axis. CDCs do not follow the positive-/7 
to negative- /7 mapping of the Fermi-Hubbard Hamilto- 
nian. CDCs decrease as /7 > is increased, since the 
system is closer to a Mott product state. In the fourth 
row of Fig. [6] we show the pairing correlations, which are 
of interest only in the negative- /7 region where pairing 
naturally occurs due to the negative interaction. 

In the fifth row of Fig. [6l we show the antiferromag- 
netic structure factor. This measure is the spatial fourier 
transform of the spin-spin correlation function evaluated 
at g = TT, and is motivated by the fact that a canonical 
transformation maps the half- filled positive- Hubbard 
model in the strong coupling limit U/t ^ 1 onto a spin- 
1/2 antiferromagnetic Heisenberg model with J = At'^/U 
in an external magnetic field h = 2/i 



H 



(18) 



This measure jumps sharply at half filling where the an- 
tiferromagnetic order is expected to exist, grows pos- 
itively and reaches a maximum for some value of /7, 
and then decreases to for large U. The symmetry 
about /i is preserved because the operators involved are 
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FIG. 6: (color online) Quantum Measures for the Fermi- Hubbard Hamiltonian in Particle-Hole Symmetric Form. Density (first 
row), generalized entanglement or Q-measure (second row), charge- density wave correlation (third row), pairing correlation 
(fourth row), antiferromagnetic structure factor (fifth row), and fidelity susceptibility (sixth row). The panels show increasing 
numbers of sites from left to right: 3, 5, 9, 10, and 20. 
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which vanish for on-site singlet 



5« - 
pairs. 

In the sixth row of Fig. [6] we show the fidelity suscep- 
tibility. Due to our use of the grand canonical ensemble, 
only different total number states have an overlap signif- 
icantly different from unity for small systems. For larger 



systems, fidelity susceptibility becomes less dominated 
by differences in number states as the number of jumps 
in number approaches the resolution in (/i/t, /7/t), as was 
also the case for the Bose-Hubbard Hamiltonian in Fig.O 
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FIG. 7: (color online) Comparison of Canonical Ensemble 
(CE) and Grand Canonical Ensemble (CCE). Shown is the 
Q-measure (generalized entanglement) for the Bose-Hubbard 
Hamiltonian and L — 20 sites. The CE is restricted to curves 
of fixed number difference. To obtain the images in Fig. [21 we 
interpolate the CE by taking the nearest fixed line and thus 
"filling in" the phase diagram. The GCE shows that this 
method obtains the correct result, thus overcoming the finite 
differencing error inherent to the CE: ji ^ E{N + 1) — E{N). 
Small qualitative differences remain between the two ensem- 
bles, as can be observed by comparing to the corresponding 
panel (second row, fourth column) in Fig. (2] 



IV. KEY NUMERICAL COMPARISONS 

A. Grand Canonical vs. Canonical Ensembles 

The canonical ensemble is much preferable for large 
systems because it is more efficient. However, in or- 
der to obtain a chemical potential one must approxi- 
mate the derivative in /i = {dE / dN)\y . Our volume 
is automatically fixed by performing the derivative for 
a given system size L. The derivative is then evaluated 
by finite differences, either a simple forward difference 
of /i ^ E{N ^ 1) — E{N)^ or a higher order finite dif- 
ference. At any level of finite difference approximation, 
for smaller systems there is only a very low resolution 
for the chemical potential axis. This leads to a lack of 
information about quantum measures. This issue is high- 
lighted in Fig. [71 by showing what the actual data is for 
the canonical ensemble (right panel) as compared to the 
grand canonical ensemble (left panel) for 20 sites. In the 
canonical calculation, each horizontal curve represents an 
allowed value of N. The grand canonical ensemble al- 
lows us to motivate the interpolation of the canonical 
ensemble result to go past the finite- difference result and 
obtain a continuous plot. Our interpolation algorithm 
is to take the value from the ne arest curve according to 
the standard distance measure vvA^^^^^^^^J^^^^i^^2)^- 
Figure [8] shows the actual canonical ensemble data used 
to produce Fig. [2l 



The initial state chosen for our imaginary time propa- 
gation has the same projection onto every Fock state in 
the Hilbert space allowed by our numerical truncation. 
In our grand canonical calculations, this includes states 
which do not have the same total number, A^. However, 
for small systems, number states are sufficiently different 
in energy that imaginary time propagation converges on 
a state with a definite TV. On the border between num- 
ber states one requires a very high x to deal with the 
near degeneracy of states with different N. This leads 
to the regular pattern of bright spots seen in the left 
panel of Fig. [3 The sole effect of higher x is to remove 
bright spots from the entropy and Q-measure. The den- 
sity of such low-entropy spikes increases with increasing 
L because a higher x is needed for a larger system - more 
sites means the possibility for more spatial entanglement. 
This is another reason why canonical ensemble calcula- 
tions are preferred for bosons. We have verified that all 
panels in Fig. [2] match the grand canonical picture, with 
the exception of the bright spots, which do not occur in 
canonical ensemble calculations. 

Although we do not illustrate it in Fig. [7J the variance 
in the total number, (7V^) — (A^)^, where iV = ^ • hf\ is 
zero to 5 or 6 digits of precision. This is in fact a measure 
of convergence for our grand canonical simulations; we 
require all Schmidt values (eigenvalues of the reduced 
density matrix) be converged to 10 digits in imaginary 
time propagation, and obtain this total number variance, 
having started with an initial state with weights on all 
total numbers allowed by the cut-off (i, i.e., from atoms 
to (d — 1) X L atoms. 



B. Gutzwiller Mean Field vs. Entangled 
Time-Evolving Block Decimation 

Since mean field methods are usually the most 
tractable approaches to solution of many body problems, 
identifying their domain of validity is important. For in- 
stance, for the three-dimensional Bose-Hubbard Hamilto- 
nian one can obtain the phase diagram via the Hubbard- 
Stratonovich transformation [35[, which is a mean field 
method, or alternately, via the Gutzwiller mean field ap- 
proximation, as discussed in Sec. H On the other hand, 
the Gutzwiller method fails for the one-dimensional Bose- 
Hubbard Hamiltonian [g^. It is natural to ask, how far 
beyond mean field does one need to go in order to obtain 
correct results for essential quantum measures? 

To address this question, we repeat the calculations of 
Sec, nil Bl for the Fermi-Hubbard Hamiltonian for 10 sites 
only, starting with the mean field calculation of x = 1 
and then doubling x successively. Figure [9] shows the 
results: each column corresponds to x = 1,2,4,8 while 
the rows correspond to the same quantum measures as 
Fig. [6l Note that x = 16 can be found in the third 
column of Fig. [6l so we do not repeat it here. We make 
several remarks concerning this figure. The density is a 
first moment and only requires lower values of x, as can 
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FIG. 8: Quantum Measures for the Bose-Hubbard Hamltonian in the Canonical Ensemble. This figure shows the actual data 
underlying Fig. [21 for smaller systems, the canonical ensemble provides very poor resolution in chemical potential, as can be 
observed in the first few columns. As the system size increases, the resolution in chemical potential goes up and the Mott lobe 
tip appears. 



be seen in the first row. For x = 1 in the first column, the 
entanglement shown in the second row is zero to 16 digits 
of accuracy, reflecting our use of double precision floating 
point arithmetic in our implementation of TEBD. The 
entanglement does not converge until x = 16, as can be 
determined from a comparison of Fig. [9] and Fig. [6l 

In the third row, for large positive U/t we obtain the 
correct result for charge density correlations at half- filling 
even for very low x because the Mott phase is nearly a 
product state. The pairing phase for negative /7/t, on the 
other hand, requires higher x, as can be seen in the fourth 
row. This can be understood from the fact that superflu- 
ids and superconductors are highly entangled in position 
space; similarly, in the Bose Hubbard Hamiltonian one 
must go to high x in order to resolve the superfluid region 
far away from the Mott lobes. By "spatial entanglement" 
we refer to exactly the Schmidt number used in TEBD. 
The fifth row shows anti- ferromagnetic order. 



Finally, in the sixth row we observe that the main 
features of the fidelity susceptibility appear already for 
X = 2, because fidelity measures in the grand canonical 
ensemble are mainly sensitive to total number states in 
our system; correct total number states are obtained at 
very low x- We observe that particle- hole symmetry is 
maintained even at the mean field level of x = 1, while 
the positive- to negative- /7 mapping at half- filling and 
zero chemical potential is only good for at least X = 2. 

Although Fig. [9] shows that x = 16 is sufficient to ob- 
tain by-eye converged values of all quantum measures 
considered in Fig. [6l a more careful quantitative study is 
desirable. In Fig. [10] we show this study for both fermions 
and bosons for 10 sites. Both the average error (left pan- 
els) and maximum error (right panels) over all pixels in 
{li/t^U/t) for fermions and {ii/t^t/U) for bosons are eval- 
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FIG. 9: (color online) Increasing Corrections to Mean Field Gutzwiller Ansatz. The same quantities as calculated in Fig. [6] for 
the Fermi-Hubbard Hamiltonian, but for 10 sites and in increasing values of x from left to right: 1,2,4,8. X = 1 corresponds 
to the Gutzwiller mean field approximation. The speckle structure observed in the vacuum region and for entanglement in the 
X = 1 case is simply numerical noise associated with the calculation of small numbers. 



uated. The error for any given pixel p is defined as 



Sp{M) 



i[M^(2^-) + M^(2^-i)] 



(19) 



our case taken to be 1, 2, 3, 4, 5. That is, we evaluate 
measure M at x of 16, 8, 4, 2 and compare it to the eval- 
uation of the same measure at the previous smaller value 



where M{2^) is the quantum measure of interest evalu- 
ated at Schmidt number x — 2-^ j is an integer, in 
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of X' The maximum and average errors are given by 

(M) = max£:p(M) , (20) 

„M) . (21, 

For fermions we calculate these two errors for density, Q- 
measure, and pairing in the grand canonical ensemble for 
10 sites, corresponding to Fig. [9l For brevity, we do not 
show an analogous figure to Fig. [9] for bosons, but calcu- 
late maximum and average error for L = 10 sites and the 
grand canonical ensemble shown in Fig. [3 in the third 
column. As can be seen from the figures, the average er- 
ror is about 1 % for X = 16 for bosons, and the maximum 
error is about 10%. For fermions it somewhat larger for 
the same x- The maximum error appears large, but upon 
a detailed investigation we find that this error is localized 
in the vacuum or band insulator regions; this means that 
the max error is dominated by imaginary time propaga- 
tion effects. Other regions have an error corresponding 
to the average error, including the important half- filling 
Mott region for positive U. In general, local measures 
like density or entanglement converge faster than non- 
local measures like quantum depletion, as can be seen 
in the figure. Also, the bosonic measures converge more 
quickly than do the fermionic measures. 

V. CONCLUSIONS 

We explored finite size effects in the Bose- and Fermi- 
Hubbard Hamiltonians with TEBD and imaginary time 
propagation, studying a range of quantum measures in- 
cluding moments, correlations, entanglement, and fi- 
delity. We found that entanglement [iol , followed 
the same pattern as more traditional measures for quan- 
tum phase transitions, such as depletion for the Bose- 
Hubbard Hamiltonian and the superfluid-Mott-insulator 
transition, and pairing, charge-density, and antiferromag- 
netic structure factor for the Fermi-Hubbard Hamilto- 
nian in the Mott transition; the exception was the tip 
of the Mott lobe, which was best described by the de- 
pletion, in terms of converging to the thermodynamic 
result more quickly as a function of system size. Fidelity 
susceptibility provides a useful alternative measure for 
characterizing the phase transition, and is a non-local 
measure. 

We compared results computed within the grand 
canonical ensemble and canonical ensembles: the chem- 
ical potential can be approximated in the canonical en- 
semble by a finite difference, a method only useful for 
large N. However, by first performing a grand canoni- 
cal calculation and then using it to motivate an interpo- 
lation of the canonical ensemble we could subsequently 
use the more efficient canonical calculations to obtain 
more highly converged results. The canonical ensemble 
is also more useful for fidelity susceptibility calculations. 
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FIG. 10: (color online) Quantitative Error Analysis for in- 
creasing entanglement (x) '^^ TEBD. Average error (left col- 
umn) and maximum error (right column) for the Fermi- 
Hubbard Hamiltonian (top row) and Bose-Hubbard Hamil- 
tonian (bottom row) as a function of increasing powers of x- 
For fermions, we show the quantum measures of density (red 
squares), Q-measure or generalized entanglement (blue trian- 
gles), and pairing correlations (green circles); for bosons, we 
show the error in the Q-measure (red squares) and quantum 
depletion (blue triangles) in the grand canonical ensemble. 
Curves are a guide to the eye only. 



which are otherwise dominated by changes in the con- 
served quantity of total number. 

We performed quantitative convergence studies and 
showed that particle-hole symmetry is maintained in 
the particle-hole symmetric form of the Fermi-Hubbard 
Hamiltonian for all values of x^ but the positive- to 
negative- mapping at half- filling and zero chemical po- 
tential is correct only for at least x = 2, meaning it is 
not satisfied in the mean field approximation. We found 
that for bosons the tip of the Mott lobe is wide-open for 
smaller system sizes of less than about 20 sites, so that 
exploring quantum phase transitions by changing chem- 
ical potential finds sharp effects while changing hopping 
find a very broad crossover. 

We acknowledge useful discussions with Ippei Dan- 
shita, Ludwig Mathey, Ryan Mishmash, Alejandro Mu- 
ramatsu, and Marcos Rigol. This work was supported 
by the National Science Foundation under Grant PHY- 
0547845 as part of the NSF CAREER program (LDC, 
MLW) and by the Aspen Center for Physics (LDC), 
by the Summer Undergraduate Research Fellow (SURF) 
program at the National Institute of Standards and Tech- 
nology (DCS, RCB), and by the NSF under Physics Fron- 
tiers Center grant PHY-0822671 (RCB, CWC). 



13 



S. Sachdev, Quantum Phase Transitions (Cambridge 

University Press, New York, 1999). [24 

M. A. Caprio, P. Cejnar, and F. lachello, Ann. Phys. [25 

323, 1106 (2008). [26^ 

K. Kim, M.-S. Chang, R. Islam, S. Korenblit, L.-M. [27 

Duan, and C. Monroe, arXiv:0905.0225 (2009). [28 
W. Hofstetter et al, Phys. Rev. Lett. 89, 220407 (2002). 

D. Jaksch and P. Zoher, Annals of Physics (2004). [29 
A. Lewenstein, M. Sanpera et a/.. Adv. Phys. 56, 243 

(2007). [3o; 

I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. [31 
80, 885 (2008). 

T. D. Kiihner, S. R. White, and H. Monien, Phys. Rev. [32 
B 61, 12474 (2000). [33^ 
G. G. Batrouni et a/., Phys. Rev. Lett. 89, 117203 (2002). 
N. Gemelke, X. Zhang, C.-L. Hung, and C. Chin, [34 
arXiv:0904.1532 (2009). 

U. Schneider et a/.. Science 322, 1520 (2008). [35; 
R. Jordens et a/.. Nature (London) 455, 204 (2008). 
M. Rigol, G. G. Batrouni, V. G. Rousseau, and R. T. [36 
Scalettar, Phys. Rev. A 79, 053605 (2009). 

M. Rigol and A. Muramatsu, Phys. Rev. A 69, 053612 [37 

(2004) . 

S. Trotzky et aL, arXiv:0905.4882 (2009). [38 
T. D. Kiihner and H. Monien, Phys. Rev. B 58, R14741 [39^ 
(1998). [40 
K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 

(2005) . 

A. Fubini et aL, e-print cond-mat/0605602 (2006). [41 

E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 
(1968). [42' 
S.-J. Gu, S.-S. Deng, Y.-Q. Li, and H.-Q. Lin, Phys. Rev. 
Lett. 93, 086402 (2004). [43 
L. C. Venuti et aL, Phys. Rev. B 78, 115410 (2008). [44^ 
L Titvinidze, M. Snoek, and W. Hofstetter, Phys. Rev. [45 
B 79, 144506 (2009). [46^ 

[23] A. J. Daley, C. Kollath, U. Schollwock, and G. Vidal, J. 



Stat. Mech. - Theor. Exp. 2004, P04005 (2004). 
| http:/ / physics.mines.edu / downloads / software / tebd/] 
O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956) . 
R. T. Scalettar et a/., Phys. Rev. Lett. 62, 1407 (1989). 

C. N. Varney et aL, arXiv:0903.2519 (2009). 

A. Kitaev and J. Preskih, Phys. Rev. Lett. 96, 110404 
(2006). 

D. A. Meyer and N. R. Wallach, Journal of Mathematical 
Physics 43, 4273 (2002). 

G. K. Brennen, Quant. Inf. Comp. 3, 619 (2003). 

H. Barnum, E. Knill, G. Ortiz, and L. Viola, Phys. Rev. 
A 68, 032308 (2003). 

H. Barnum et a/., Phys. Rev. Lett. 92, 107902 (2004). 
P. Buonsante and A. Vezzani, Phys. Rev. Lett. 98, 

110601 (2007). 

M. Rigol, B. S. Shastry, and S. Haas, larXiv: 0906. 5347] 
(2009). 

M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. 
Fisher, Phys. Rev. B 40, 546 (1989). 
H. Monien, M. Linn, and N. Elstner, Phys. Rev. A 58, 
R3395 (1998). 

L. D. Carr and M. K. Oberthaler, Phys. Rev. Lett., under 
review (2009). 

U. Schollwock, Rev. Mo d. Phys. 77, 259 (2005). 
[http://geco.mines.edu/ 



M. A. Nielsen and I. L. Chuang, Quantum Computation 
and Quantum Information (Cambridge University Press, 
Cambridge, UK, 2000). 

A. Osterloh, L. Amico, G. Falci, and R. Fazio, Nature 
416, 608 (2002). 

G. O. S. Deng and L. Viola, e-print 

'http://arxiv.org/abs/0802.3941 (2008). 

R. Somma et a/., Phys. Rev. A 70, 042311 (2004). 

P. W. Anderson, Phys. Rev. Lett. 18, 1049 (1967). 

E. H. Lieb and F. Y. Wu, Physica A 321, . 

T.-C. Wei et a/., Phys. Rev. A 71, 060305 (2005). 



